Entrega I

Author

Eva López Corazón - DNI: 71313437-M

Instrucciones (leer antes de empezar)

  • Modifica de la cabecera del documento .qmd tus datos personales (nombre y DNI). IMPORTANTE: no modifiques nada más de la cabecera.

  • Asegúrate, ANTES de seguir editando el documento, que el archivo .qmd se renderiza correctamente y se genera el .html correspondiente en tu carpeta local de tu ordenador (pulsa el botón Render o guarda el documento con Render on Save activado).

  • Los chunks (cajas de código) creados están o vacíos o incompletos, de ahí que la mayoría tengan la opción #| eval: false. Una vez que edites lo que consideres, debes ir cambiando cada chunck a #| eval: true (o quitarlo directamente) para que se ejecuten.

  • Recuerda que puedes ejecutar chunk a chunk con el botón play o ejecutar todos los chunk hasta uno dado (con el botón a la izquierda del anterior).

  • IMPORTANTE: solo se podrá subir un archivo .html al campus, no se evaluarán entregas sin el .html correctamente generado. Recuerda: el código es una herramienta, no el fin en esta asignatura. Se evaluará especialmente las interpretaciones y conclusiones que detalles. Un ejercicio con el código perfecto pero sin ningún tipo de razonamiento, interpretación o conclusión no superará el 4 sobre 10 de nota.

Paquetes necesarios

Necesitaremos los siguientes paquetes (haz play en el chunk para que se carguen):

rm(list = ls()) # Borramos variables de environment

# descomentar si es la primera vez (y requieren instalación)
# install.packages("tidyverse")
# install.packages("performance")
# install.packages("olsrr")
# install.packages("corrr")
# install.packages("corrplot")
# install.packages("skimr")
library(tidyverse)
library(rsample)
library(performance)
library(olsrr)
library(corrr)
library(corrplot)
library(skimr)

Caso práctico: análisis de datos de cáncer

El archivo de datos a usar lo cargaremos desde el csv cancer.csv.

# no cambies código
datos <- read_csv(file = "./cancer.csv")
datos
# A tibble: 3,047 × 13
   deathRate medIncome popEst2015 povertyPercent studyPerCap MedianAge
       <dbl>     <dbl>      <dbl>          <dbl>       <dbl>     <dbl>
 1      165.     61898     260131           11.2       500.       39.3
 2      161.     48127      43269           18.6        23.1      33  
 3      175.     49348      21026           14.6        47.6      45  
 4      195.     44243      75882           17.1       343.       42.8
 5      144.     49955      10321           12.5         0        48.3
 6      176      52313      61023           15.6       180.       45.4
 7      176.     37782      41516           23.2         0        42.6
 8      184.     40189      20848           17.8         0        51.7
 9      190.     42579      13088           22.3         0        49.3
10      178.     60397     843954           13.1       428.       35.8
# ℹ 3,037 more rows
# ℹ 7 more variables: AvgHouseholdSize <dbl>, PercentMarried <dbl>,
#   PctHS25_Over <dbl>, PctUnemployed16_Over <dbl>, PctPrivateCoverage <dbl>,
#   BirthRate <dbl>, region <chr>

Los datos representan diferentes características socioeconómicas de distintas regiones, extraídas de the American Community Survey (census.gov), clinicaltrials.gov, y cancer.gov.

¿Nuestro objetivo? Ser capaces de predecir de manera lineal y MULTIVARIANTE la variable deathRate, nuestra variable objetivo, que representa la mortalidad media de cancer, por cada 100 000 habitantes. El resto de variables son:

  • medianIncome: mediana de los ingresos de la región.
  • popEst2015: población de la región
  • povertyPercent: porcentaje de población en situación de pobreza.
  • studyPerCap: ensayos clínicos relacionados por el cáncer realizados por cada 100 000 habitantes.
  • MedianAge: edad mediana.
  • region: nombre de la región.
  • AvgHouseholdSize: tamaño medio de los hogares.
  • PercentMarried: porcentaje de habitantes casados.
  • PctHS25_Over: porcentaje de residentes por encima de los 25 años con (máximo) título de bachillerato.
  • PctUnemployed16_Over: porcentaje de residentes de más de 16 años en paro.
  • PctPrivateCoverage: porcentaje de residentes con seguro de salud privado.
  • BirthRate: tasa de natalidad.

IMPORTANTE: siempre que puedas, relaciona los valores numéricos de la salida de R con su fórmula.

Ejercicio 1 (1.25 puntos)

Chequea que no hay problemas de rango/codificación: en caso de que alguna variable tenga dichos problemas, debes sustituir los valores fueran del rango razonable por su media o mediana (sin esos datos fuera del rango), justificando por qué usas la media o por qué usas la mediana, ejecutando el código que consideres (si te atascas, sigue haciendo con el dataset original)

datos |> str()
spc_tbl_ [3,047 × 13] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
 $ deathRate           : num [1:3047] 165 161 175 195 144 ...
 $ medIncome           : num [1:3047] 61898 48127 49348 44243 49955 ...
 $ popEst2015          : num [1:3047] 260131 43269 21026 75882 10321 ...
 $ povertyPercent      : num [1:3047] 11.2 18.6 14.6 17.1 12.5 15.6 23.2 17.8 22.3 13.1 ...
 $ studyPerCap         : num [1:3047] 499.7 23.1 47.6 342.6 0 ...
 $ MedianAge           : num [1:3047] 39.3 33 45 42.8 48.3 45.4 42.6 51.7 49.3 35.8 ...
 $ AvgHouseholdSize    : num [1:3047] 2.54 2.34 2.62 2.52 2.34 2.58 2.42 2.24 2.38 2.65 ...
 $ PercentMarried      : num [1:3047] 52.5 44.5 54.2 52.7 57.8 50.4 54.1 52.7 55.9 50 ...
 $ PctHS25_Over        : num [1:3047] 23.2 26 29 31.6 33.4 30.4 29.8 31.6 32.2 28.8 ...
 $ PctUnemployed16_Over: num [1:3047] 8 7.8 7 12.1 4.8 12.9 8.9 8.9 10.3 9.2 ...
 $ PctPrivateCoverage  : num [1:3047] 75.1 70.2 63.7 58.4 61.6 60 49.5 55.8 55.5 69.9 ...
 $ BirthRate           : num [1:3047] 6.12 4.33 3.73 4.6 6.8 ...
 $ region              : chr [1:3047] "West" "West" "West" "West" ...
 - attr(*, "spec")=
  .. cols(
  ..   deathRate = col_double(),
  ..   medIncome = col_double(),
  ..   popEst2015 = col_double(),
  ..   povertyPercent = col_double(),
  ..   studyPerCap = col_double(),
  ..   MedianAge = col_double(),
  ..   AvgHouseholdSize = col_double(),
  ..   PercentMarried = col_double(),
  ..   PctHS25_Over = col_double(),
  ..   PctUnemployed16_Over = col_double(),
  ..   PctPrivateCoverage = col_double(),
  ..   BirthRate = col_double(),
  ..   region = col_character()
  .. )
 - attr(*, "problems")=<externalptr> 
datos |> skim()
Data summary
Name datos
Number of rows 3047
Number of columns 13
_______________________
Column type frequency:
character 1
numeric 12
________________________
Group variables None

Variable type: character

skim_variable n_missing complete_rate min max empty n_unique whitespace
region 0 1 4 9 0 4 0

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
deathRate 0 1 178.66 27.75 59.70 161.20 178.10 195.20 362.80 ▁▇▆▁▁
medIncome 0 1 47063.28 12040.09 22640.00 38882.50 45207.00 52492.00 125635.00 ▇▇▁▁▁
popEst2015 0 1 102637.37 329059.22 827.00 11684.00 26643.00 68671.00 10170292.00 ▇▁▁▁▁
povertyPercent 0 1 16.88 6.41 3.20 12.15 15.90 20.40 47.40 ▃▇▃▁▁
studyPerCap 0 1 155.40 529.63 0.00 0.00 0.00 83.65 9762.31 ▇▁▁▁▁
MedianAge 0 1 45.27 45.30 22.30 37.70 41.00 44.00 624.00 ▇▁▁▁▁
AvgHouseholdSize 0 1 2.48 0.43 0.02 2.37 2.50 2.63 3.97 ▁▁▃▇▁
PercentMarried 0 1 51.77 6.90 23.10 47.75 52.40 56.40 72.50 ▁▂▇▇▁
PctHS25_Over 0 1 34.80 7.03 7.50 30.40 35.30 39.65 54.80 ▁▂▇▇▁
PctUnemployed16_Over 0 1 7.85 3.45 0.40 5.50 7.60 9.70 29.40 ▅▇▁▁▁
PctPrivateCoverage 0 1 64.35 10.65 22.30 57.20 65.10 72.10 92.30 ▁▂▇▇▂
BirthRate 0 1 5.64 1.99 0.00 4.52 5.38 6.49 21.33 ▂▇▁▁▁ 
#region variable cualitativa
#variable edad media se pasa
unique(datos$MedianAge) |> sort()
  [1]  22.3  23.2  23.3  23.5  23.9  24.2  24.4  24.6  24.7  25.0  25.2  25.9
 [13]  26.1  26.2  26.3  26.4  26.5  26.6  26.7  26.8  27.0  27.2  27.4  27.5
 [25]  27.6  27.8  27.9  28.0  28.1  28.3  28.4  28.5  28.6  28.7  28.8  28.9
 [37]  29.1  29.2  29.3  29.4  29.5  29.7  29.8  29.9  30.0  30.1  30.2  30.3
 [49]  30.4  30.5  30.6  30.7  30.8  30.9  31.0  31.1  31.2  31.3  31.4  31.5
 [61]  31.6  31.7  31.8  31.9  32.0  32.1  32.2  32.3  32.4  32.5  32.6  32.7
 [73]  32.8  32.9  33.0  33.1  33.2  33.3  33.4  33.5  33.6  33.7  33.8  33.9
 [85]  34.0  34.1  34.2  34.3  34.4  34.5  34.6  34.7  34.8  34.9  35.0  35.1
 [97]  35.2  35.3  35.4  35.5  35.6  35.7  35.8  35.9  36.0  36.1  36.2  36.3
[109]  36.4  36.5  36.6  36.7  36.8  36.9  37.0  37.1  37.2  37.3  37.4  37.5
[121]  37.6  37.7  37.8  37.9  38.0  38.1  38.2  38.3  38.4  38.5  38.6  38.7
[133]  38.8  38.9  39.0  39.1  39.2  39.3  39.4  39.5  39.6  39.7  39.8  39.9
[145]  40.0  40.1  40.2  40.3  40.4  40.5  40.6  40.7  40.8  40.9  41.0  41.1
[157]  41.2  41.3  41.4  41.5  41.6  41.7  41.8  41.9  42.0  42.1  42.2  42.3
[169]  42.4  42.5  42.6  42.7  42.8  42.9  43.0  43.1  43.2  43.3  43.4  43.5
[181]  43.6  43.7  43.8  43.9  44.0  44.1  44.2  44.3  44.4  44.5  44.6  44.7
[193]  44.8  44.9  45.0  45.1  45.2  45.3  45.4  45.5  45.6  45.7  45.8  45.9
[205]  46.0  46.1  46.2  46.3  46.4  46.5  46.6  46.7  46.8  46.9  47.0  47.1
[217]  47.2  47.3  47.4  47.5  47.6  47.7  47.8  47.9  48.0  48.1  48.2  48.3
[229]  48.4  48.5  48.6  48.7  48.8  48.9  49.0  49.1  49.2  49.3  49.4  49.5
[241]  49.6  49.7  49.8  49.9  50.0  50.1  50.2  50.3  50.4  50.5  50.6  50.7
[253]  50.8  50.9  51.0  51.1  51.2  51.3  51.4  51.5  51.7  51.8  51.9  52.0
[265]  52.1  52.2  52.3  52.4  52.5  52.6  52.7  52.8  52.9  53.1  53.2  53.3
[277]  53.4  53.5  53.6  53.7  54.0  54.1  54.4  54.5  54.6  54.7  54.8  55.1
[289]  55.4  55.6  56.3  56.4  56.5  56.6  57.1  57.3  59.0  65.3 349.2 406.8
[301] 412.8 414.0 424.8 429.6 430.8 458.4 469.2 470.4 481.2 496.8 498.0 500.4
[313] 501.6 502.8 508.8 511.2 519.6 523.2 525.6 535.2 536.4 546.0 579.6 619.2
[325] 624.0
datos_preproc <- 
  datos |> 
  mutate(MedianAge = ifelse(MedianAge > 130, median(MedianAge), MedianAge))
unique(datos_preproc$MedianAge) |> sort()
  [1] 22.3 23.2 23.3 23.5 23.9 24.2 24.4 24.6 24.7 25.0 25.2 25.9 26.1 26.2 26.3
 [16] 26.4 26.5 26.6 26.7 26.8 27.0 27.2 27.4 27.5 27.6 27.8 27.9 28.0 28.1 28.3
 [31] 28.4 28.5 28.6 28.7 28.8 28.9 29.1 29.2 29.3 29.4 29.5 29.7 29.8 29.9 30.0
 [46] 30.1 30.2 30.3 30.4 30.5 30.6 30.7 30.8 30.9 31.0 31.1 31.2 31.3 31.4 31.5
 [61] 31.6 31.7 31.8 31.9 32.0 32.1 32.2 32.3 32.4 32.5 32.6 32.7 32.8 32.9 33.0
 [76] 33.1 33.2 33.3 33.4 33.5 33.6 33.7 33.8 33.9 34.0 34.1 34.2 34.3 34.4 34.5
 [91] 34.6 34.7 34.8 34.9 35.0 35.1 35.2 35.3 35.4 35.5 35.6 35.7 35.8 35.9 36.0
[106] 36.1 36.2 36.3 36.4 36.5 36.6 36.7 36.8 36.9 37.0 37.1 37.2 37.3 37.4 37.5
[121] 37.6 37.7 37.8 37.9 38.0 38.1 38.2 38.3 38.4 38.5 38.6 38.7 38.8 38.9 39.0
[136] 39.1 39.2 39.3 39.4 39.5 39.6 39.7 39.8 39.9 40.0 40.1 40.2 40.3 40.4 40.5
[151] 40.6 40.7 40.8 40.9 41.0 41.1 41.2 41.3 41.4 41.5 41.6 41.7 41.8 41.9 42.0
[166] 42.1 42.2 42.3 42.4 42.5 42.6 42.7 42.8 42.9 43.0 43.1 43.2 43.3 43.4 43.5
[181] 43.6 43.7 43.8 43.9 44.0 44.1 44.2 44.3 44.4 44.5 44.6 44.7 44.8 44.9 45.0
[196] 45.1 45.2 45.3 45.4 45.5 45.6 45.7 45.8 45.9 46.0 46.1 46.2 46.3 46.4 46.5
[211] 46.6 46.7 46.8 46.9 47.0 47.1 47.2 47.3 47.4 47.5 47.6 47.7 47.8 47.9 48.0
[226] 48.1 48.2 48.3 48.4 48.5 48.6 48.7 48.8 48.9 49.0 49.1 49.2 49.3 49.4 49.5
[241] 49.6 49.7 49.8 49.9 50.0 50.1 50.2 50.3 50.4 50.5 50.6 50.7 50.8 50.9 51.0
[256] 51.1 51.2 51.3 51.4 51.5 51.7 51.8 51.9 52.0 52.1 52.2 52.3 52.4 52.5 52.6
[271] 52.7 52.8 52.9 53.1 53.2 53.3 53.4 53.5 53.6 53.7 54.0 54.1 54.4 54.5 54.6
[286] 54.7 54.8 55.1 55.4 55.6 56.3 56.4 56.5 56.6 57.1 57.3 59.0 65.3

- `deathRate`: entre 0 e infinito.

- `medianIncome`: entre 0 e infinito.

- `popEst2015`: entre 0 y 8.000.000.000 (habitantes del mundo xd)

- `povertyPercent`: entre 0 y 100 (al ser porcentaje)

  • `studyPerCap`: entre 0 e infinito

  • `MedianAge`: entre 0 y 130 años (por poner un tope)

  • `AvgHouseholdSize`: entre 0 e infinito (si existe alguna finca enorme)

- `PercentMarried`: entre 0 y 100 (al ser porcentaje)

- `PctHS25_Over`: entre 0 y 100 (al ser porcentaje)

  • `PctUnemployed16_Over`: entre 0 y 100 (al ser porcentaje)

- `ctPrivateCoverage`: entre 0 y 100 (al ser porcentaje)

- `BirthRate`: entre 0 e infinito.

a la hora de sustituir usa lo mediana porque es menos sensible a valores extremos como un 500 y pico que había

Ejercicio 2 (0.75 puntos)

En caso de encontrar alguna variable cualitativa, procesa los datos como consideras para incluir su información de manera que pueda ser usada en la futura regresión. Recuerda: si tenemos una variable cualitativa que toma k valores, debemos crear k-1 nuevas variables numéricas (de un tipo muy concreto) en su lugar.

# region cualitativa
unique(datos$region)
[1] "West"      "South"     "Midwest"   "Northeast"
# hay 4 valores posibles: creamos 3 variables nuevas binarias
datos_preproc <-
  datos_preproc |> 
  mutate('West' = ifelse(region == "West", 1, 0),
            'South' = ifelse(region == "South", 1, 0),
         'Midwest' = ifelse(region == "Midwest", 1, 0)) |> 
  select(-region)

region es cualitativa con 4 posibles calores: creamso 3 nuevas variables y quitamos region

Ejercicio 3 (1.5 puntos)

Ejecuta el código que consideres para realizar una selección inicial de variables de manera que se mitiguen los efectos adversos de la colinealidad. Comenta todo lo que consideres sobre. Tras ello separa las observaciones en dos datasets distintos: un dataset train y un segundo dataset test. Para ello realiza un muestreo aleatorio de manera que train contenga el 80% de los datos y test el 20% restante. IMPORTANTE: no cambies la semilla.

#correlacion con deathrate
#datos_preproc |> correlate() |> focus(deathRate)
#datos_preproc |> cor() |> corrplot(method='ellipse')
#quitamos predictoras muy pococ correladas para simplificarnos los calculos
datos_2<- datos_preproc |> dplyr::select(-c(studyPerCap,MedianAge,AvgHouseholdSize))


set.seed(12345)
split<-initial_split(datos_2, prop=0.8 )
train<-training(split)
test<-testing(split)

por tomar una medida quitamos las variables poco correladas con la objetivo deathrate, en este caso las menores al |0.1| aunque tambien podriamos considerar quitar popEst2015

Ejercicio 4 (1 punto)

No usaremos el dataset test hasta el final para evaluar las predicciones (con un dataset que el modelo no conoce). Con el dataset train ejecuta el ajuste de regresión lineal multivariante saturado (habiendo hecho lo pedido en los ejercicios anteriores) y comenta de manera MUY RESUMIDA la salida (SOLO lo que puedas interpretar hasta ahora, y al grano, que se te acaba el tiempo)

ajuste_saturado<-lm(train, formula=deathRate~.)
ajuste_saturado |> summary()

Call:
lm(formula = deathRate ~ ., data = train)

Residuals:
     Min       1Q   Median       3Q      Max 
-108.344  -11.959    0.627   11.555  158.574 

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)    
(Intercept)           1.243e+02  1.210e+01  10.280  < 2e-16 ***
medIncome            -3.574e-05  7.566e-05  -0.472 0.636724    
popEst2015           -9.133e-07  1.350e-06  -0.676 0.498794    
povertyPercent        6.144e-01  1.812e-01   3.391 0.000708 ***
PercentMarried       -3.148e-01  9.774e-02  -3.221 0.001296 ** 
PctHS25_Over          1.267e+00  8.404e-02  15.079  < 2e-16 ***
PctUnemployed16_Over  1.551e+00  1.863e-01   8.324  < 2e-16 ***
PctPrivateCoverage    6.610e-02  8.496e-02   0.778 0.436627    
BirthRate            -5.372e-01  2.338e-01  -2.297 0.021683 *  
West                 -8.133e+00  2.195e+00  -3.704 0.000217 ***
South                 1.038e+01  1.953e+00   5.317 1.15e-07 ***
Midwest               3.136e+00  1.956e+00   1.603 0.109030    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 21.94 on 2425 degrees of freedom
Multiple R-squared:  0.3765,    Adjusted R-squared:  0.3737 
F-statistic: 133.1 on 11 and 2425 DF,  p-value: < 2.2e-16
# comprobamos colinealidad entre predictoras con VIF
check_collinearity(ajuste_saturado)
# Check for Multicollinearity

Low Correlation

                 Term  VIF   VIF 95% CI Increased SE Tolerance Tolerance 95% CI
            medIncome 4.33 [4.04, 4.64]         2.08      0.23     [0.22, 0.25]
           popEst2015 1.19 [1.15, 1.26]         1.09      0.84     [0.80, 0.87]
       PercentMarried 2.28 [2.15, 2.43]         1.51      0.44     [0.41, 0.47]
         PctHS25_Over 1.75 [1.66, 1.86]         1.32      0.57     [0.54, 0.60]
 PctUnemployed16_Over 2.00 [1.89, 2.12]         1.41      0.50     [0.47, 0.53]
   PctPrivateCoverage 4.17 [3.89, 4.47]         2.04      0.24     [0.22, 0.26]
            BirthRate 1.07 [1.04, 1.14]         1.04      0.93     [0.88, 0.96]
                 West 3.07 [2.87, 3.28]         1.75      0.33     [0.31, 0.35]
                South 4.78 [4.46, 5.13]         2.19      0.21     [0.19, 0.22]
              Midwest 4.29 [4.01, 4.61]         2.07      0.23     [0.22, 0.25]

Moderate Correlation

           Term  VIF   VIF 95% CI Increased SE Tolerance Tolerance 95% CI
 povertyPercent 6.78 [6.30, 7.30]         2.60      0.15     [0.14, 0.16]
train2<- train|> select(-povertyPercent)

ajuste_saturado2<-lm(train2, formula=deathRate~.)
ajuste_saturado2 |> summary()

Call:
lm(formula = deathRate ~ ., data = train2)

Residuals:
     Min       1Q   Median       3Q      Max 
-103.288  -11.893    0.604   11.863  158.564 

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)    
(Intercept)           1.569e+02  7.375e+00  21.275  < 2e-16 ***
medIncome            -1.690e-04  6.480e-05  -2.608 0.009175 ** 
popEst2015           -1.060e-06  1.352e-06  -0.784 0.433138    
PercentMarried       -4.777e-01  8.531e-02  -5.599 2.40e-08 ***
PctHS25_Over          1.228e+00  8.342e-02  14.721  < 2e-16 ***
PctUnemployed16_Over  1.669e+00  1.834e-01   9.096  < 2e-16 ***
PctPrivateCoverage   -5.346e-02  7.747e-02  -0.690 0.490157    
BirthRate            -5.073e-01  2.342e-01  -2.167 0.030357 *  
West                 -8.090e+00  2.200e+00  -3.677 0.000241 ***
South                 1.117e+01  1.943e+00   5.745 1.03e-08 ***
Midwest               3.588e+00  1.956e+00   1.835 0.066667 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 21.99 on 2426 degrees of freedom
Multiple R-squared:  0.3736,    Adjusted R-squared:  0.371 
F-statistic: 144.7 on 10 and 2426 DF,  p-value: < 2.2e-16
check_collinearity(ajuste_saturado2)
# Check for Multicollinearity

Low Correlation

                 Term  VIF   VIF 95% CI Increased SE Tolerance Tolerance 95% CI
            medIncome 3.16 [2.96, 3.38]         1.78      0.32     [0.30, 0.34]
           popEst2015 1.19 [1.14, 1.25]         1.09      0.84     [0.80, 0.87]
       PercentMarried 1.73 [1.64, 1.83]         1.32      0.58     [0.55, 0.61]
         PctHS25_Over 1.72 [1.63, 1.82]         1.31      0.58     [0.55, 0.61]
 PctUnemployed16_Over 1.93 [1.82, 2.05]         1.39      0.52     [0.49, 0.55]
   PctPrivateCoverage 3.45 [3.23, 3.69]         1.86      0.29     [0.27, 0.31]
            BirthRate 1.07 [1.04, 1.13]         1.03      0.93     [0.88, 0.96]
                 West 3.07 [2.87, 3.28]         1.75      0.33     [0.31, 0.35]
                South 4.71 [4.39, 5.06]         2.17      0.21     [0.20, 0.23]
              Midwest 4.27 [3.99, 4.59]         2.07      0.23     [0.22, 0.25]

consideramos un VIF mayor a 5 demasiado por lo que eliminamos povertyPercent, recomprobamos y parece que las variables predictoras no estan significativamente correladas entre si

Ejercicio 5 (1.75 puntos)

Realiza una selección de modelos en base a los criterios BIC y AIC. Recuerda que para evitar incompatibilidad con otros paquetes, no debes hAcer library(MASS) sino MASS::… Comenta e interpreta todo lo que sepas de las salidas generadas. ¿Cuál penaliza más? ¿Por qué? ¿Qué ventajas tiene la selección de BIC? Interpreta los coeficientes de ambos modelos y comenta las diferencias (en caso de que existiesen)

#BIC
MASS::stepAIC(ajuste_saturado2, direction='both',k=log(nrow(train2)))
Start:  AIC=15137.43
deathRate ~ medIncome + popEst2015 + PercentMarried + PctHS25_Over + 
    PctUnemployed16_Over + PctPrivateCoverage + BirthRate + West + 
    South + Midwest

                       Df Sum of Sq     RSS   AIC
- PctPrivateCoverage    1       230 1172938 15130
- popEst2015            1       297 1173005 15130
- Midwest               1      1627 1174335 15133
- BirthRate             1      2269 1174977 15134
- medIncome             1      3287 1175995 15136
<none>                              1172708 15137
- West                  1      6535 1179243 15143
- PercentMarried        1     15155 1187863 15161
- South                 1     15955 1188663 15163
- PctUnemployed16_Over  1     39994 1212702 15211
- PctHS25_Over          1    104760 1277467 15338

Step:  AIC=15130.11
deathRate ~ medIncome + popEst2015 + PercentMarried + PctHS25_Over + 
    PctUnemployed16_Over + BirthRate + West + South + Midwest

                       Df Sum of Sq     RSS   AIC
- popEst2015            1       251 1173189 15123
- Midwest               1      1537 1174475 15126
- BirthRate             1      2125 1175064 15127
<none>                              1172938 15130
- West                  1      6310 1179248 15135
- medIncome             1      6658 1179596 15136
+ PctPrivateCoverage    1       230 1172708 15137
- PercentMarried        1     15328 1188266 15154
- South                 1     16986 1189924 15157
- PctUnemployed16_Over  1     48377 1221315 15221
- PctHS25_Over          1    104776 1277714 15331

Step:  AIC=15122.84
deathRate ~ medIncome + PercentMarried + PctHS25_Over + PctUnemployed16_Over + 
    BirthRate + West + South + Midwest

                       Df Sum of Sq     RSS   AIC
- Midwest               1      1623 1174812 15118
- BirthRate             1      2108 1175297 15119
<none>                              1173189 15123
- West                  1      6201 1179390 15128
- medIncome             1      7272 1180460 15130
+ popEst2015            1       251 1172938 15130
+ PctPrivateCoverage    1       184 1173005 15130
- PercentMarried        1     15080 1188269 15146
- South                 1     17441 1190630 15151
- PctUnemployed16_Over  1     48126 1221315 15213
- PctHS25_Over          1    109109 1282298 15332

Step:  AIC=15118.41
deathRate ~ medIncome + PercentMarried + PctHS25_Over + PctUnemployed16_Over + 
    BirthRate + West + South

                       Df Sum of Sq     RSS   AIC
- BirthRate             1      1742 1176554 15114
<none>                              1174812 15118
+ Midwest               1      1623 1173189 15123
+ popEst2015            1       337 1174475 15126
+ PctPrivateCoverage    1        98 1174714 15126
- medIncome             1      9521 1184333 15130
- PercentMarried        1     13708 1188520 15139
- West                  1     25366 1200178 15163
- South                 1     31422 1206234 15175
- PctUnemployed16_Over  1     47178 1221990 15207
- PctHS25_Over          1    107638 1282450 15324

Step:  AIC=15114.22
deathRate ~ medIncome + PercentMarried + PctHS25_Over + PctUnemployed16_Over + 
    West + South

                       Df Sum of Sq     RSS   AIC
<none>                              1176554 15114
+ BirthRate             1      1742 1174812 15118
+ Midwest               1      1257 1175297 15119
+ popEst2015            1       308 1176246 15121
+ PctPrivateCoverage    1        26 1176528 15122
- medIncome             1      8811 1185365 15125
- PercentMarried        1     15444 1191998 15138
- West                  1     25515 1202069 15159
- South                 1     32512 1209067 15173
- PctUnemployed16_Over  1     47741 1224296 15203
- PctHS25_Over          1    109486 1286040 15323

Call:
lm(formula = deathRate ~ medIncome + PercentMarried + PctHS25_Over + 
    PctUnemployed16_Over + West + South, data = train2)

Coefficients:
         (Intercept)             medIncome        PercentMarried  
           1.541e+02            -2.141e-04            -4.624e-01  
        PctHS25_Over  PctUnemployed16_Over                  West  
           1.233e+00             1.694e+00            -1.072e+01  
               South  
           8.706e+00  
# salida ultima iteracion:lm(formula = deathRate ~ medIncome + PercentMarried + PctHS25_Over + PctUnemployed16_Over + West + South, data = train2)
ajuste_BIC <-lm(formula = deathRate ~ medIncome + PercentMarried + PctHS25_Over + 
    PctUnemployed16_Over + West + South, 
    data = train2)
ajuste_BIC |> summary()

Call:
lm(formula = deathRate ~ medIncome + PercentMarried + PctHS25_Over + 
    PctUnemployed16_Over + West + South, data = train2)

Residuals:
     Min       1Q   Median       3Q      Max 
-105.425  -12.091    0.501   11.766  159.447 

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)    
(Intercept)           1.541e+02  5.944e+00  25.933  < 2e-16 ***
medIncome            -2.141e-04  5.020e-05  -4.266 2.07e-05 ***
PercentMarried       -4.624e-01  8.187e-02  -5.648 1.82e-08 ***
PctHS25_Over          1.233e+00  8.202e-02  15.038  < 2e-16 ***
PctUnemployed16_Over  1.694e+00  1.706e-01   9.930  < 2e-16 ***
West                 -1.072e+01  1.476e+00  -7.259 5.21e-13 ***
South                 8.706e+00  1.062e+00   8.194 4.02e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 22 on 2430 degrees of freedom
Multiple R-squared:  0.3715,    Adjusted R-squared:   0.37 
F-statistic: 239.4 on 6 and 2430 DF,  p-value: < 2.2e-16
#AIC
MASS::stepAIC(ajuste_saturado2, direction='both',k=2)
Start:  AIC=15073.65
deathRate ~ medIncome + popEst2015 + PercentMarried + PctHS25_Over + 
    PctUnemployed16_Over + PctPrivateCoverage + BirthRate + West + 
    South + Midwest

                       Df Sum of Sq     RSS   AIC
- PctPrivateCoverage    1       230 1172938 15072
- popEst2015            1       297 1173005 15072
<none>                              1172708 15074
- Midwest               1      1627 1174335 15075
- BirthRate             1      2269 1174977 15076
- medIncome             1      3287 1175995 15078
- West                  1      6535 1179243 15085
- PercentMarried        1     15155 1187863 15103
- South                 1     15955 1188663 15105
- PctUnemployed16_Over  1     39994 1212702 15153
- PctHS25_Over          1    104760 1277467 15280

Step:  AIC=15072.13
deathRate ~ medIncome + popEst2015 + PercentMarried + PctHS25_Over + 
    PctUnemployed16_Over + BirthRate + West + South + Midwest

                       Df Sum of Sq     RSS   AIC
- popEst2015            1       251 1173189 15071
<none>                              1172938 15072
- Midwest               1      1537 1174475 15073
+ PctPrivateCoverage    1       230 1172708 15074
- BirthRate             1      2125 1175064 15074
- West                  1      6310 1179248 15083
- medIncome             1      6658 1179596 15084
- PercentMarried        1     15328 1188266 15102
- South                 1     16986 1189924 15105
- PctUnemployed16_Over  1     48377 1221315 15169
- PctHS25_Over          1    104776 1277714 15279

Step:  AIC=15070.65
deathRate ~ medIncome + PercentMarried + PctHS25_Over + PctUnemployed16_Over + 
    BirthRate + West + South + Midwest

                       Df Sum of Sq     RSS   AIC
<none>                              1173189 15071
- Midwest               1      1623 1174812 15072
+ popEst2015            1       251 1172938 15072
+ PctPrivateCoverage    1       184 1173005 15072
- BirthRate             1      2108 1175297 15073
- West                  1      6201 1179390 15082
- medIncome             1      7272 1180460 15084
- PercentMarried        1     15080 1188269 15100
- South                 1     17441 1190630 15105
- PctUnemployed16_Over  1     48126 1221315 15167
- PctHS25_Over          1    109109 1282298 15285

Call:
lm(formula = deathRate ~ medIncome + PercentMarried + PctHS25_Over + 
    PctUnemployed16_Over + BirthRate + West + South + Midwest, 
    data = train2)

Coefficients:
         (Intercept)             medIncome        PercentMarried  
           1.536e+02            -2.012e-04            -4.713e-01  
        PctHS25_Over  PctUnemployed16_Over             BirthRate  
           1.238e+00             1.705e+00            -4.851e-01  
                West                 South               Midwest  
          -7.770e+00             1.149e+01             3.564e+00  
# salida ultima iteracion:lm(formula = deathRate ~ medIncome + PercentMarried + PctHS25_Over + PctUnemployed16_Over + BirthRate + West + South + Midwest, data = train2)
ajuste_AIC <- lm(formula = deathRate ~ medIncome + PercentMarried + PctHS25_Over + PctUnemployed16_Over + BirthRate + West + South + Midwest, 
    data = train2)
ajuste_AIC |> summary()

Call:
lm(formula = deathRate ~ medIncome + PercentMarried + PctHS25_Over + 
    PctUnemployed16_Over + BirthRate + West + South + Midwest, 
    data = train2)

Residuals:
     Min       1Q   Median       3Q      Max 
-103.450  -11.926    0.524   11.870  158.525 

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)    
(Intercept)           1.536e+02  6.226e+00  24.673  < 2e-16 ***
medIncome            -2.012e-04  5.187e-05  -3.879 0.000108 ***
PercentMarried       -4.713e-01  8.436e-02  -5.587 2.58e-08 ***
PctHS25_Over          1.238e+00  8.236e-02  15.027  < 2e-16 ***
PctUnemployed16_Over  1.705e+00  1.709e-01   9.980  < 2e-16 ***
BirthRate            -4.851e-01  2.323e-01  -2.089 0.036829 *  
West                 -7.770e+00  2.169e+00  -3.582 0.000347 ***
South                 1.149e+01  1.912e+00   6.008 2.16e-09 ***
Midwest               3.564e+00  1.945e+00   1.833 0.066953 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 21.98 on 2428 degrees of freedom
Multiple R-squared:  0.3733,    Adjusted R-squared:  0.3712 
F-statistic: 180.8 on 8 and 2428 DF,  p-value: < 2.2e-16

BIC:

deathRate = 1.541e+02+-2.141e-04*medIncome + -4.624e-01*PercentMarried + 1.233e+00*PctHS25_Over +1.694e+00* PctUnemployed16_Over + -1.072e+01*West + 8.706e+00*South

AIC:

deathRate =           1.536e+02 +
medIncome          *  -2.012e-04  +
PercentMarried      * -4.713e-01  +
PctHS25_Over         * 1.238e+00  +
PctUnemployed16_Over *  1.705e+00  +
BirthRate           * -4.851e-01  + 
West                * -7.770e+00  +
South               *  1.149e+01  +
Midwest             *  3.564e+00  +

Como vemos enl BIC penaliza más, quita dos predictoras mas que AIC, las predictoras birthrate y midwest , que no son buenas: cambia precision por senicillez.

Ejercicio 6 (1.5 puntos)

Compara la calidad de los 3 modelos (saturado, BIC y AIC) en train. Quédate con uno (justifica) y comprueba si se cumplen las hipotesis necesarias (fase de diagnosis SOLO en ese elegido). ¿Se cumplen las hipótesis? Argumenta porqué, más allá de lo que valga luego la bondad de ajuste, es importante que se cumplan.

compare_performance(ajuste_saturado2, ajuste_AIC,ajuste_BIC)
# Comparison of Model Performance Indices

Name             | Model |   AIC (weights) |  AICc (weights) |   BIC (weights) |    R2 | R2 (adj.) |   RMSE |  Sigma
--------------------------------------------------------------------------------------------------------------------
ajuste_saturado2 |    lm | 21991.6 (0.154) | 21991.7 (0.151) | 22061.1 (<.001) | 0.374 |     0.371 | 21.936 | 21.986
ajuste_AIC       |    lm | 21988.6 (0.690) | 21988.6 (0.691) | 22046.5 (0.013) | 0.373 |     0.371 | 21.941 | 21.982
ajuste_BIC       |    lm | 21991.5 (0.156) | 21991.6 (0.158) | 22037.9 (0.987) | 0.372 |     0.370 | 21.972 | 22.004

Obviamente el AIC ofrece menor AIC y el BIC menor BIC y el saturado presenta el mayor AIC y BIC.

Como es de esperar tambien el saturado ofrece mayor R2 (infromación eplicada) que AIC y BIC y AIC mayor R2 que BIC.

El BIC penaliza más, con menos variables pero el AIC ofrece mas informacion (apenas) a mas parametros aunque estos sean malos y por lo general es preferible elegir un modelo mas simple, especialmente la diferencia de informacion explicada es tan pequeña (0.002 ó 0.001 en R2 adj).

–>BIC es el mejor en este caso (sencillo a similar precisión)

Es necesario realizar inferencia estadística delos estimadores de los parámetros para comprobar que no solo funcionan en esta muestra, si no en otra y se pueda aplicar a toda la población

  1. test linealidad -> OK

    check_model(ajuste_BIC)

    chisq.test(ajuste_BIC$residuals,ajuste_BIC$fitted.values, simulate.p.value = TRUE, B=250)
    
    Pearson's Chi-squared test with simulated p-value (based on 250
    replicates)

data:  ajuste_BIC$residuals and ajuste_BIC$fitted.values
X-squared = 5936532, df = NA, p-value = 0.003984

  2. test homocedasticidad -> NO

    necesitamos que la varianza del error sea finita y constante, tal que no varíe según aumenta o disminuye la X.

    performance::check_heteroscedasticity(ajuste_BIC)
    Warning: Heteroscedasticity (non-constant error variance) detected (p < .001).

  3. test normalidad de erroes -> NO

    ε∼N(0,σ2r)

performance::check_normality(ajuste_BIC)
Warning: Non-normality of residuals detected (p < .001).

  1. independencia de errores -> OK

    performance::check_autocorrelation(ajuste_BIC)
    OK: Residuals appear to be independent and not autocorrelated (p = 0.622).

Ejercicio 7 (0.75 puntos)

Asume que se cumplen las hipótesis (en caso de que no se cumpliesen) e interpreta todo lo que sepas sobre la segunda tabla de la salida de la regresión (la de inferencia de los parámetros). Elimina variables si su efecto no fuese significativo.

ajuste_saturado2 |> summary()

Call:
lm(formula = deathRate ~ ., data = train2)

Residuals:
     Min       1Q   Median       3Q      Max 
-103.288  -11.893    0.604   11.863  158.564 

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)    
(Intercept)           1.569e+02  7.375e+00  21.275  < 2e-16 ***
medIncome            -1.690e-04  6.480e-05  -2.608 0.009175 ** 
popEst2015           -1.060e-06  1.352e-06  -0.784 0.433138    
PercentMarried       -4.777e-01  8.531e-02  -5.599 2.40e-08 ***
PctHS25_Over          1.228e+00  8.342e-02  14.721  < 2e-16 ***
PctUnemployed16_Over  1.669e+00  1.834e-01   9.096  < 2e-16 ***
PctPrivateCoverage   -5.346e-02  7.747e-02  -0.690 0.490157    
BirthRate            -5.073e-01  2.342e-01  -2.167 0.030357 *  
West                 -8.090e+00  2.200e+00  -3.677 0.000241 ***
South                 1.117e+01  1.943e+00   5.745 1.03e-08 ***
Midwest               3.588e+00  1.956e+00   1.835 0.066667 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 21.99 on 2426 degrees of freedom
Multiple R-squared:  0.3736,    Adjusted R-squared:  0.371 
F-statistic: 144.7 on 10 and 2426 DF,  p-value: < 2.2e-16

Ahora que se ‘cumplen las hipótesis’ podemos interpretar el p-valor, R2 y otros estadisticos

Ejercicio 8 (1.5 puntos)

Por último, utiliza el dataset test para, aplicando los 3 modelos entrenados (saturado, AIC y BIC), predecir las observaciones de test (con cada uno de los modelos. Con las observaciones de test construye un nuevo dataset de TRES columnas: la variable objetivo, la variable predictora, la predicción y el modelo usado (saturado, BIC, AIC). Calcula el \(R^2\) en el dataset de test para cada modelo. ¿Cuál funciona mejor en test? ¿Qué % de información explica cada modelo? ¿Cuánto vale su varianza residual?

pred_test_S <-
  tibble("objetivo" = test$deathRate,
         "prediccion" = predict(ajuste_saturado2, newdata = test),
         "split" = "test")
pred_test_A <-
  tibble("objetivo" = test$deathRate,
         "prediccion" = predict(ajuste_AIC, newdata = test),
         "split" = "test")
pred_test_B <-
  tibble("objetivo" = test$deathRate,
         "prediccion" = predict(ajuste_BIC, newdata = test),
         "split" = "test")

pred_train_S <-
  tibble("objetivo" = train$deathRate,
         "prediccion" = ajuste_saturado2$fitted.values,
         "split" = "train")
pred_train_A <-
  tibble("objetivo" = train$deathRate,
         "prediccion" = ajuste_AIC$fitted.values,
         "split" = "train")
pred_train_B <-
  tibble("objetivo" = train$deathRate,
         "prediccion" = ajuste_BIC$fitted.values,
         "split" = "train")
prediccion <- bind_rows(pred_test_S, pred_test_A, pred_test_B, pred_train_S, pred_train_A, pred_train_B)

prediccion
# A tibble: 9,141 × 3
   objetivo prediccion split
      <dbl>      <dbl> <chr>
 1     144.       155. test 
 2     237.       206. test 
 3     207.       191. test 
 4     210.       208. test 
 5     156.       183. test 
 6     140.       193. test 
 7     169.       199. test 
 8     158.       159. test 
 9     176.       185. test 
10     187.       173. test 
# ℹ 9,131 more rows